library(tidyverse)
library(gWQS)
library(ggplot2)
library(cowplot)
library(knitr)
library(kableExtra)
library(gtsummary)
library(haven)
library(skimr)
library(corrplot)
library(bkmr)
library(forestplot)
library(mice)
library(rms)
library(tableone)
library(ggplot2)
library(ggcorrplot)
library(survminer)
library(sjmisc)
library(meta)
library(dplyr)
library(ggplot2)
library(survminer)
library(survival)
library(forestmodel)

# setwd('G:/BaiduNetdiskDownload/nhanes')
# 
# # 年龄 种族 BMI 
# # 年龄-RIDAGEYR; 性别-RIAGENDR; 种族-RIDRETH3; 教育程度-DMDEDUC2; 贫困程度-INDFMPIR;
# # 血清性类固醇激素 
# # 尿砷代谢物
# # 基础的年龄性别种族
# demo_h <- read_xpt("2013-2014/Demographics/demo_h.xpt")
# demo_i <- read_xpt("2015-2016/Demographics/demo_i.xpt")
# demo_all <- dplyr::bind_rows(list(demo_h,demo_i))
# # 可替宁
# cot_h <- read_xpt("2013-2014/Laboratory/cot_h.xpt")
# cot_i <- read_xpt("2015-2016/Laboratory/cot_i.xpt")
# cot_all <- dplyr::bind_rows(list(cot_h,cot_i)) 
# # BMI 
# BMI_h <- read_xpt("2013-2014/Examination/bmx_h.xpt")
# BMI_i <- read_xpt("2015-2016/Examination/bmx_i.xpt")
# bmi_all <- dplyr::bind_rows(list(BMI_h,BMI_i))
# 
# #Urinary Arsenic 尿砷
# # Laboratory	UTAS	urxuas	Urinary arsenic, total (ug/L)	尿砷总量 (ug/L)	Urinary arsenic, total (ug/L)
# 
# utas_h <- read_xpt("2013-2014/Laboratory/utas_h.xpt")
# utas_i <- read_xpt("2015-2016/Laboratory/utas_i.xpt")
# utas_all <- dplyr::bind_rows(list(utas_h,utas_i))
# 
# #Urinary Arsenic  metabolites 尿砷代谢物  Urinary creatinine (mg/dL) 也在这里面
# # Arsenic, total
# # 检测了三种尿砷代谢物  亚砷酸 (Arsenous acid UAS	urxuas3)、DMA (二甲基胂酸 UAS urxudma)、MMA(一甲基胂酸 UAS	urxumma)
# uas_h <- read_xpt("2013-2014/Laboratory/uas_h.xpt")
# uas_i <- read_xpt("2015-2016/Laboratory/uas_i.xpt")
# uas_all <- dplyr::bind_rows(list(uas_h,uas_i))
# # 三种血清性类固醇激素
# # Total testosterone (TT ) 总睾酮 TST	lbxtst,  estradiol (E2) 雌二醇 Laboratory	TST	lbxest , sex hormone-binding globulin   性激素结合球蛋白(SHBG) Laboratory	TST	lbxshbg
# tst_h <- read_xpt("2013-2014/Laboratory/tst_h.xpt")
# tst_i <- read_xpt("2015-2016/Laboratory/tst_i.xpt")
# tst_all <- dplyr::bind_rows(list(tst_h,tst_i))
# 
# output <- plyr::join_all(list(demo_all, cot_all,uas_all,utas_all,tst_all,bmi_all),by='SEQN',type='full')
# table(output$URXUAS3)
# 
# # 根据图筛选6-19岁
# output619 <- output[output$RIDAGEYR >=6 &output$RIDAGEYR<=19, ]
# # 正确 5451
# dim(output619) #5451
# 
# # 排除不含有尿砷的人 3826 我这边3817
# output619Arsenous <- output619[which(!is.na(output619$URXUAS)), ]
# # 排除缺失性类固醇激素
# 
# output619Arsenouslbx <- output619Arsenous[which(!is.na(output619Arsenous$LBXTST)), ]
# output619Arsenouslbx <- output619Arsenouslbx[which(!is.na(output619Arsenouslbx$LBXEST)), ]
# output619Arsenouslbx <- output619Arsenouslbx[which(!is.na(output619Arsenouslbx$LBXSHBG)), ]
# 
# # 排除 BMI,  BMI cate, 血清可替宁, 收入,
# output619Arsenouslbx <- output619Arsenouslbx[which(!is.na(output619Arsenouslbx$BMXBMI)), ]
# output619Arsenouslbx <- output619Arsenouslbx[which(!is.na(output619Arsenouslbx$BMDBMIC)), ]
# output619Arsenouslbx <- output619Arsenouslbx[which(!is.na(output619Arsenouslbx$LBXCOT)), ]
# output619Arsenouslbx <- output619Arsenouslbx[which(!is.na(output619Arsenouslbx$INDFMPIR)), ]

dim(output619Arsenouslbx)
# 数据写入本地进行减少内存使用
# write.csv(output619Arsenouslbx,file = 'F:/Rproject/r-language/paper2/output619Arsenouslbx.csv')


#tab3 表格测试森林图 做不了 转移出去
table(table1datademoAllParse$Sex)
# 分出四个亚组分区数据来
tab3dataMaleChildren <-subset(table1datademoAllParse,Sex=='Male'&childgroup=='Children')
tab3dataMaleAdolescents <-subset(table1datademoAllParse,Sex=='Male'&childgroup=='Adolescents')
tab3dataFeMaleChildren <-subset(table1datademoAllParse,Sex=='Female'&childgroup=='Children')
tab3dataFeMaleChildren <-subset(table1datademoAllParse,Sex=='Female'&childgroup=='Adolescents')

tab3dataMaleChildrenm1 <- glm(Arsenictotal ~ TT+Age+race+BMI+PIR+Cotinineexposurestatus,
                              data =tab3dataMaleChildren)




tab3dataMaleChildrenm1reg <-tbl_regression(tab3dataMaleChildrenm1,exponentiate = T,include =c(TT))


show_header_names(tab3dataMaleChildrenm1reg)

tab3dataMaletab3dataMaleAdolescents <- glm(Arsenictotal ~ TT+Age+race+BMI+PIR+Cotinineexposurestatus,
                                           data =tab3dataMaleAdolescents)

tab3dataMaleChildrenm1Adolescents <-tbl_regression(tab3dataMaletab3dataMaleAdolescents,exponentiate = T,include =c(TT))

mmm1 <-tbl_stack(tbls = list(tab3dataMaleChildrenm1reg,tab3dataMaleChildrenm1Adolescents),group_header = c('male','male'))
F1 <-forest_model(tab3dataMaletab3dataMaleAdolescents,covariates  = 'TT',)

pretty_lung <- lung %>%
  transmute(time,
            status,
            Age = age,
            Sex = factor(sex, labels = c("Male", "Female")),
            ECOG = factor(lung$ph.ecog),
            `Meal Cal` = meal.cal)
write.csv(table1datademoAllParse,file = 'aaa.csv')

# forest_model(tab3dataMaleChildrenm1,covariates  = 'TT')
# 2023年7月19日11:07:59 查阅了大量文章发现不好做 要么是每个模型生成数据 保存在excel 用excel进行画出森林图